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PROJECT  SUMMARY 


The  hypersonic,  viscous,  chemically  reacting  flow  around  a  vehicle  in  the 
upper  layers  of  the  atmosphere  influences  vehicle  performance,  thermal 
environment,  and  sensor  visibility.  Under  the  Phase  I  effort,  a  unique  hybrid 
Navier- Stokes/Monte  Carlo  method  has  been  successfully  developed  for  tne 
simulation  of  said  hypersonic  flows.  This  method  combines  a  continuum  description 
of  the  gas  mixture  with  a  statistical  (or  "particle")  description  of  the  species 
transport  and  the  chemical  reactions ,  thereby  allowing  easy  use  of  available 
chemical  kinetics  data.  As  such,  the  hybrid  approach  is  capable  of  exploiting 
both  the  efficiency  of  a  continuum  approach  for  the  fluid  dynamic  calculations  and 
the  benefits  in  physical  modeling  of  the  Monte  Carlo  approach  for  chemical 
reactions.  The  viability  of  this  hybrid  approach  for  practical  hypersonic  flow 
calculations  has  been  demonstrated  in  the  Phase  I  effort  via  the  calculation  of 
the  chemically  reacting  flow  field  over  an  axisymmetric  cone.  A  proposed  Phase  II 
effort  would  continue  the  development  of  this  method  by  improving  the  numerical 
efficiency  of  the  particle  transport  scheme,  and  by  incorporating  physical  models 
for  catalytic  walls,  recombination  reactions,  and  ionization.  In  addition, 
radiation  would  be  included  via  a  radiation  transport  model.  For  the  class  of 
problems  of  interest  here,  the  resulting  procedure  will  be  able  to  represent  the 
complex  chemical  and  thermal  processes  occurring  around  a  hypersonic  vehicle  in 
much  more  detail  than  is  currently  possible  with  either  the  continuum  or  the 
direct  simulation  Monte  Carlo  (DSMC)  approach.  The  proposed  Phase  II  program 
would  also  modify  the  computer  code  for  use  in  a  workstation  environment  via  the 
development  of  suitable  input  and  output  packages,  thus  producing  a  valuable 
design  tool  for  the  examination  of  the  thermal  environment  and  the  sensor 
visibility  of  hypersonic  vehicles. 
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1 .  INTRODUCTION 


Ground  or  air  launched  interceptors  can  achieve  high  velocities  (>  3  kft/sec) 
within  the  atmosphere  either  during  the  boost  or  intercept  phase  of  their  mission. 
Space  launched  interceptors  may  also  find  it  advantageous  to  be  able  to  penetrate 
the  atmosphere  to  enlarge  their  domain  of  effectiveness.  The  high  velocities 
necessary  to  achieve  successful  interception  result  in  high  Mach  numbers  (>  10) 
within  the  atmosphere  (<  300  kft) .  The  resulting  frictional  heating  of  the 
boundary  layer  and  the  collisional  heating  due  to  the  bow  shock  create  gas 
temperatures  that  are  sufficiently  high  (>  6000  K)  that  significant  energy  is 
absorbed  in  disassociating  the  air  and  significant  energy  transport  can  occur  by 
radiation.  This  radiation  is  complicated  by  the  fact  that  the  gas  is  partly 
disassociated.  The  disassociation  and  radiation  causes  two  specific  concerns. 

The  first  is  the  subsequent  additional  thermal  loading  over  and  above  the  normal 
convective  load  on  the  thermal  protection  system  due  to  recombination  and 
radiative  transfer  effects  which,  if  known,  could  be  accounted  for  in  the 
design.  The  second  concern  is  potentially  more  important  since  it  bears  on  the 
vehicle  functionality  and  relates  to  the  vehicle  signature  and  the  ability  of 
sensors  to  see  source  radiant  energy  through  this  complex,  vehicle- induced 
chemical  and  radiative  environment. 

Insofar  as  this  source  radiant  energy  is  concerned,  it  is  not  simply  the 
integrated  flux  at  the  sensor  that  is  required,  but  the  distribution  across  the 
spectrum  is  of  interest,  for  purposes  of  sensor  design  and  source  identification. 
The  transmission  of  this  source  radiant  energy  through  the  vehicle  environment 
unfortunately  depends  on  the  gas  composition  and  temperature,  since  the  absorption 
and  emissivity  are  functions  of  both.  As  a  result,  certain  frequencies  suffer 
much  more  absorption  than  others,  particularly  at  the  lower  temperatures  and  at 
frequencies  around  4000A  (.4^m),  the  near  infrared,  but  the  extent  is  a  very 
complex  issue  requiring  very  detailed  information  on  the  local  gas  composition  and 
temperature . 

Potential  additional  complication  in  an  already  complex  problem  arises  with 
active  cooling  of  the  vehicle  boundary  layer  or  ablation.  Such  coolant  would 
reduce  or  change  the  radiation  from  the  boundary  layer,  but  also  would  absorb 
radiation  from  the  shock.  Note  also  that  the  boundary  layer,  the  coolant  flow  and 
the  shock  structure  are  very  strongly  interacting  in  hypersonic  flow  and  thus  a 
coupled  flow,  chemistry  analysis  and  radiation  would  seem  appropriate.  The  degree 
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to  which  the  vehicle  thermal  environment  as  well  as  radiation  at  the  sensor  would 
change  due  to  coolant  flow  is  difficult  to  estimate  without  doing  a  very  precise 
computation. 

It  is  further  observed  that  the  shock  wave  structure,  particularly  the  bow 
structure,  is  very  important  in  setting  the  flow  and  thermal  fields  (via  the 
temperature),  the  radiation  field  and  the  degree  of  disassociation  (which  also 
influences  the  temperature),  which  in  turn  greatly  influences  the  spectral  dis¬ 
tribution  of  radiation  emissivity  and  adsorptivity .  This  shock  structure,  loca¬ 
tion,  and  strength  in  terms  of  pressure  and  temperature  rise  is  in  turn  determined 
by  the  forebody  geometry,  incidence,  flight  velocity  and  altitude.  Thus,  forebody 
geometry  at  a  particular  flight  condition  could  be  examined  and  varied  to  see  if  a 
beneficial  impact  on  signature,  thermal  environment  and  sensor  visibility  could  be 
achieved.  However,  in  order  to  make  such  an  assessment,  very  detailed,  coupled 
reacting  flow  and  radiation  field  predictions  are  necessary.  It  is  the 
development  of  a  procedure  for  efficiently  making  these  calculations  and  the 
performance  of  these  predictions  that  is  the  subject  of  the  present  Phase  I  and 
the  proposed  Phase  II  effort.  Some  of  the  complexities  of  the  flow  field  around 
this  type  of  vehicle  are  illustrated  in  the  volume  edited  by  Nelson  (Ref.  1). 


2 .  BACKGROUND 

The  major  task  that  has  to  be  performed  is  the  development  of  a  hybrid  viscous 
flow  field  Computational  Fluid  Dynamics  (CFD)  prediction  code  which  would  be 
suitable  for  efficiently  determining  the  vehicle  flow  and  thermal  environment, 
including  the  effects  of  disassociation  and  radiation,  at  high  Mach  numbers  within 
the  atmosphere.  During  the  Phase  I  effort  the  major  emphasis  has  been  placed  upon 
demonstrating  the  potential  of  a  unique  hybrid  Navier- Stokes/Monte  Carlo  method 
with  particular  emphasis  upon  disassociation.  The  proposed  Phase  II  effort  would 
continue  the  Phase  I  disassociation  work  and  also  incorporate  a  radiation 
capability  within  this  unique  procedure.  Before  proceeding,  it  is  useful  to 
consider  possible  approaches  to  solving  the  high  altitude,  high  Mach  number 
external  flow  problem.  These  approaches  include  usual  Navier-Stokes  procedures, 
usual  Monte -Carlo  procedures  and  the  unique  procedure  demonstrated  in  the  Phase 
I  effort.  The  following  discussion  indicates  the  advantages  of  the  Phase  I 
approach  for  the  class  of  problems  of  interest  here. 
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In  regard  to  candidate  codes,  most  existing  three-dimensional  Navier-Stokes 
codes  could  be  modified  to  include  radiative  transport  by  adding  the  appropriate 
terra  in  the  energy  equation  and  incorporating  a  radiative  transport  equation.  The 
choice  of  the  most  suitable  Navier-Stokes  code  is  not  primarily  determined  by 
radiative  effects.  The  situation  in  regard  to  chemical  reaction  is  different.  In 
recent  years,  there  has  been  a  considerable  increase  in  the  number  of  potentially 
applicable  reacting  flow  codes,  such  as  the  axisymmetric  codes  of  Candler  and 
MacCormack  (Ref.  2)  and  Green  (Ref.  3),  for  instance.  While  codes  such  as  those 
of  Refs.  2  and  3  are  based  on  the  continuum  flow  equations,  even  conceptually  more 
general  codes  based  on  the  direct  simulation  of  molecular  dynamics  using  the 
kinetic  theory  of  gases  have  been  constructed,  for  instance  by  Moss  and  Bird  (Ref. 
4).  While  such  direct  simulation  codes  are  of  great  theoretical  interest,  they 
are  currently  not  practical  in  terms  of  flow  dimensionality  and  computer  run  time 
for  the  near  continuum  and  continuum  flow  regime  and  will  require  much  development 
before  they  can  address  the  present  problem  in  any  detail.  Appropriate  codes 
based  on  the  continuum  equations  are  much  further  along  in  the  development, 
although,  for  instance,  the  codes  of  Refs.  2  and  3  are  currently  restricted  to 
axisymmetric  flow,  but  this  restriction  could  be  easily  removed. 

At  the  proposed  flight  conditions  the  flow  field  becomes  disassociated  behind 
the  bow  shock  and  the  temperature  field  reflects  the  energy  required  to 
disassociate,  and  in  some  instances  ionize,  the  air.  In  the  subsequent 
recombination  of  the  species  some  of  this  energy  is  given  up,  thus  the  time  scale 
of  the  relaxation  can  have  a  major  influence  on  the  thermal  protection  system, 
since  the  relaxation  process  governs  the  degree  and  location  of  the  heat  release. 
The  molecules  are  also  vibrationally  and  rotationally  excited  and  this  can  also 
serve  to  store  energy.  Thus  a  realistic  simulation  of  the  flow  field  requires  a 
large  number  of  distinct  species  be  tracked.  Depending  on  the  quantities  of 
interest,  chemists  have  determined  that  as  few  as  five  species  or  as  many  as 
thirty  or  more  may  be  required  (Refs.  4,  5).  To  this  is  added  one  grouping  for 
each  of  the  excited  states  that  are  required,  including  the  ionized  species.  In 
the  conventional  continuum  approach,  as  exemplified  by  schemes  such  as  those  in 
Refs.  2  and  3,  the  resulting  set  requires  one  partial  differential  equation  for 
each  species  or  state  (Ref.  6) .  Very  quickly  the  labor  of  solving  the  chemistry 
greatly  exceeds  that  of  solving  the  fluid  conservation  equations.  For  instance, 
in  two  dimensions,  in  the  typical  block  implicit  approach,  such  as  in  Ref.  2,  five 
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species  and  three  vibrationally  excited  molecules  are  followed,  resulting  in  an 
(ll/4)^-fold  increase  in  run  time  relative  to  solving  for  a  perfect  gas,  i.e., 

p  _  run  time,  reacting  _  fNg  —  1  +  Nf  1  ^ 
run  time,  perfect  gas  L  jjf  -* 

where  Nf  is  the  number  of  equations  in  the  fluids  system  (four  in  two  dimensions, 
five  in  three  dimensions),  Ng  is  the  number  of  distinct  species  or  states,  and  F 
represents  the  ratio  of  run  times  (note  -  Ref.  2  gives  F  as  proportional  to  the 
square  of  the  number  of  species,  however,  standard  textbooks  give  the  operation 
count  of  a  block  tridiagonal  system  as  going  as  the  cube  of  the  block  size). 

There  are  various  strategies  for  reducing  the  run  time  of  these  implicit 
approaches,  attempting  to  make  it  more  nearly  proportional  to  the  ratio 
(Ng— 1  +  Nf)/Nf,  rather  than  the  cube  of  this  ratio.  Depending  on  the  specific 
reactions  involved,  these  strategies  may  or  may  not  be  successful,  conditional  on 
the  stiffness  of  the  system  and  its  size.  Nonetheless,  at  present  no  clearly 
superior  general  uncoupled  alternative  has  yet  appeared. 

In  free  molecular  and  near  free  molecular  flow  it  is  computationally  feasible 
to  follow  individual  molecules  and  determine  the  appropriate  fluxes  by  averaging 
over  many  trials.  Such  molecule -by-molecule  simulation  using  the  kinetic  theory 
of  gases  contains  a  random  selection  process  which  has  become  associated  with  the 
name  'Monte  Carlo'  methodologies.  Within  the  Monte  Carlo  framework  there  are  a 
multitude  of  different  equations  and  methods.  In  gas  dynamics  a  particular 
variant  has  been  termed  the  Direct  Simulation  Monte  Carlo  (DSMC)  method.  The  DSMC 
method  is  described  in  detail  in  the  book  by  Bird  (Ref.  7) ,  and  involves  following 
groups  of  molecules  rather  than  individual  molecules.  The  various  groups  are 
interacted  in  a  probabilistic  manner.  Thus  the  approach  may  be  thought  of  as 
representing  the  gas  by  a  smaller  than  actual  number  of  real  molecules  or 
alternatively  as  a  model  gas  with  the  actual  number  of  larger  model  gas  molecules. 
Bird  (Refs.  7,  8)  has  shown  how  DSMC  can  be  extended  to  chemically  reacting  flow 
using  the  principles  of  the  kinetic  theory  of  gases  in  a  very  straightforward 
manner.  Here  the  computational  run  time  goes  only  (linearly)  as  the  increased 
number  of  molecules  required  to  statistically  describe  the  extended  system.  Bird 
also  points  out  that  the  overall  scheme  becomes  computationally  more  and  more 
expensive  as  the  Knudsen  number  becomes  smaller,  such  that  below  Kn^ . 1  the  method 
is  probably  unworkable  in  two  or  three  space  dimensions.  Here  the  problem  is  the 
requirement  that  the  typical  cell  size  be  much  smaller  than  the  mean  free  path. 
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Thus  there  is  a  dilemma.  The  well  understood  implicit  continuum  methods  for 
reacting  flow  have  a  computational  effort  that  increases  as  the  cube  of  the  number 
of  species,  but  to  obtain  reasonable  computer  run  times  from  DSMC,  where  the  run 
time  goes  linearly  with  the  number  of  molecules,  one  requires  a  very  dilute  gas 
with  a  cell  size  small  compared  to  the  mean  free  path.  The  present  effort 
suggests  an  innovative  alternative  hybrid  'two-phase  flow'  approach.  The  mean 
motion  of  the  gas  is  described  in  an  Eulerian  frame,  the  'continuum  gas  phase'  of 
the  'two  phase  flow'.  The  simulated  molecular  phase  is  described  in  a  Lagrangian 
frame,  the  'solid  particle'  phase  of  the  two  phase  flow.  In  the  usual  two  phase 
flow  manner  the  particles  are  convected  by  the  mean  flow.  However,  the  particles 
are  now  interpreted  in  the  manner  of  Direct  Simulation  Monte  Carlo  as  groups  of 
molecules.  The  various  groups  are  interacted  in  the  same  probabalistic  manner 
used  in  the  DSMC  method,  including  the  chemical  interactions.  The  mixture 
continuum  gas  phase  is  then  recomputed  using  the  species  concentrations  resulting 
from  the  particulate  phase.  Either  at  each  time  step  of  a  transient  problem  or 
upon  reaching  steady  state  as  the  asymptote  of  an  unsteady  problem,  the  gas  and 
particulate  phases  would  be  iterated  until  both  'phases'  were  in  computational 
equilibrium.  Chemical  and  thermodynamic  equilibrium  would  not  be  a  requirement, 
however,  other  than  in  the  implied  Maxwellian  form  for  the  random  molecular 
velocity  in  the  continuum  equations. 

Intuitively,  a  large  enough  chemical  system  could  always  be  chosen  such  that 
the  hybrid  approach  would  eventually  be  more  economical  than  solving  a  large 
system  of  coupled  partial  differential  equations.  However,  within  the  continuum 
approach  strategies  do  exist  for  treating  the  coupled  chemical -fluid  system  in  an 
uncoupled  manner.  From  a  computational  chemistry  viewpoint,  a  completely 
uncoupled  explicit  treatment  of  the  chemical  system  is  not  general  and  usually  is 
very  inefficient  with  the  normal  stiff  chemical  systems.  Indeed,  the  chemistry 
problem  created  the  impetus  for  the  whole  field  of  stiff  ordinary  differential 
equation  system  solvers.  Depending  on  how  uncoupled  the  system  is  physically  and 
to  what  extent  it  has  been  assumed  uncoupled  computationally,  the  degree  of 
iteration  necessary  to  reintroduce  the  system  coupling  into  the  uncoupled 
computation  may  be  very  significant.  The  current  view  is  that  sufficient 
physical  coupling  exists  for  the  case  where  the  air  chemistry  around  the  vehicle 
is  of  interest,  to  require  the  coupled  system  approach  which,  as  previously  noted, 
at  best  has  a  labor  that  goes  as  the  cube  of  the  system  size.  The  hybrid  approach 
outlined  here  does  introduce  similar  coupling  concerns,  since  the  chemistry  and 
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fluid  are  essentially  uncoupled  over  each  time  step  or  pseudo  time  step/iteration. 
Thus,  it  is  highly  desirable  to  reduce  the  labor  of  both  the  direct  molecular 
simulation  step  and  the  fluids  step  as  much  as  possible,  so  that  the  labor  of 
iteratively  coupling  the  two  systems  will  not  be  great.  Much  effort  by  a  large 
number  of  people  is  currently  underway  to  reduce  the  required  fluids  computational 
effort  for  a  multitude  of  other  reasons,  and  Scientific  Research  Associates,  Inc. 
is  actively  following  the  literature  and  participating  in  these  developments.  As 
a  result,  it  is  only  necessary  here  to  consider  the  efficiency  of  the  Monte  Carlo 
part  of  the  calculation. 

The  hybrid  approach  developed  in  the  Phase  I  investigation  is  capable  of 
exploiting  both  the  efficiency  of  a  continuum  approach  for  the  fluid  dynamic 
calculation  of  the  gas  mixture  and  the  benefits  in  physical  modeling  of  the  Monte 
Carlo  approach  for  the  chemical  reactions.  For  the  class  of  problems  of  interest 
here,  this  approach  will  be  able  to  represent  the  complex  chemical  and  thermal 
processes  more  efficiently  than  a  DSMC  approach  and  more  accurately  than  a 
conventional  Navier-Stokes  procedure. 


3.  THE  HYBRID  NAVIER- STOKES/MONTE  CARLO  METHOD 

The  hybrid  method  developed  under  the  present  effort  allows  the  analysis  of 
chemically  reacting  gas  flows  around  a  vehicle  traveling  at  hypersonic  speeds 
within  the  atmosphere.  In  this  method,  the  Navier-Stokes  equations  are  used  to 
model  the  gas  flow,  while  a  Monte  Carlo- like  technique  is  used  to  simulate  the 
chemical  reactions.  For  each  small  time  interval  Atjjj,  the  following  steps  are 
taken; 

(1)  The  Navier-Stokes  equations  are  solved  numerically  for  the  mean  flow 
field,  together  with  the  global  continuity  and  energy  equations  for  the 
gas  mixture,  assuming  "frozen"  species  concentrations. 

(2)  Computational  particles  representing  the  species  molecules  are 
transported,  using  convection  velocities  computed  in  (1)  and  taking  into 
account  molecular  diffusion  of  the  species. 
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(?'  A  representative  set  of  reactions  is  computed  among  the  particles  by 

means  of  DSMC-like  techniques.  The  post-reaction  particle  distribution 
is  used  to  update  the  species  concentrations . 

A  detailed  description  of  each  of  the  steps  of  the  hybrid  Navier- Stokes/Monte 
Carlo  method  is  presented  below. 


3.1.  Step  1,  -  The  Continuum  Phase 

In  -he  continuum  phase,  the  mean  flow  field  calculations  are  performed  via 
numerical  solution  of  the  Navier-Stokes  equations,  and  the  global  continuity  and 
energy  equations  for  the  gas  mixture.  In  vecto’'  form  this  set  of  equations  can  be 
written  as 


dp 

—  -f-  V»pU  =  0 

at 

(1) 

a  (pu) 

'  +  V*  (pUU)  =  -Vp  +  V*T 

a  t 

(2) 

^  +  V.  (pUh)  =  ^  ^  -  ^*<1  +  ^ 

at  Dt 

(3) 

where  p  is  the  density  of  the  mixture,  U  is  the  velocity,  p  is  the  pressure,  r  is 
the  stress  tensor  (molecular  and  turbulent),  h  is  the  mixture  enthalpy  (including 
the  heat  of  formation),  $  is  the  viscous  dissipation  per  unit  volume,  q  is  the 
heat  flux  vector,  and  r  is  the  radiation  term.  D/Dt  denotes  a  material 
derivative,  i.e.  D/Dt  =  d/dt  +  U*V . 

It  should  be  noted  that  the  use  of  the  Navier-Stokes  equations  with  its  usual 
stress  strain  relationship  is  consistent  with  a  Maxwellian  distribution  of  the 
molecular  random  velocities  (see  Ref.  7). 

The  above  equations  for  the  gas  mixture  have  to  be  supplemented  by  the 
equation  of  state  for  an  ideal  gas  mixture. 
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where  R  is  the  universal  gas  constant,  T  is  the  temperature,  N  is  the  number  of 
species  in  the  mixt  :e,  Xj^  is  the  mass  fraction  of  species  i,  Xj  is  the  mole 
fraction  of  species  i  and  is  the  molecular  weight  of  species  i. 

The  temperature  T  is  related  to  the  mixture  enthalpy  h  via  the  caloric 
equation  of  state 


h 


(5) 


where  hf^  is  the  formation  enthalpy  of  species  i  at  the  temperature  ,  and  where 
‘^pi  the  specific  heat  at  constant  pressure  of  species  i.  Values  of  hj^  and  Tf, 
and  empirical  relations  between  Cp  and  T  for  different  species  can  be  found  in 
Refs.  9  and  10. 

To  complete  the  set  of  equations  (1)  -  (5),  formulas  are  needed  for  the 
mixture  viscosity  /i  (which  appears  in  the  stress  tensor)  and  the  mixture  thermal 
conductivity  k  (which  appears  in  the  heat  flux  vector).  By  Wilke's  semi-empirical 
formula  (Ref.  11), 


where  is  the  viscosity  of  species  i,  and  where  <I>£j  is  given  by 


The  viscosity  of  a  single  species  can  be  determined  from  the  kinetic  theory 
expression  /'.efs.  12,  13) 
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(MT) 

^single  species  ~  266.93  10  ^  (8) 

where  fi  is  given  in  g/cra-s,  T  in  K,  and  where  o  is  a  molecular  collision  diameter 
(in  Angstrom),  and  0^  is  a  viscosity  collision  integral,  which  can  be  approximated 
by  (Ref.  14) 


fiy  "=  1.147 


■  T 

-0.145 

r  T 

.  T,  . 

+ 

+  0-5 

1  T,  J 

_2 


(9) 


Here  is  a  molecular  temperature  diameter.  Values  for  o  and  can  be  found  in 
Ref.  15. 

Finally,  the  thermal  conductivity  k  is  given  by  Eucken's  approximation  (Ref. 

13), 


N 

k  =.I 


Xj  kj 


1=1  N 


(10) 


where  kj^  is  the  thermal  conductivity  of  species  i. 


(C 


pi 


4  Mi  ' 


Mi 


(11) 


and  where  is  the  same  function  as  in  Eqs.  (6)  -  (7). 

At  M=15  and  250  kft  altitude,  the  freestream  Reynolds  number  per  foot  is 
3.5x10-’.  Flight  experiments  give  the  transition  Reynolds  number  on  cones  at  these 
Mach  numbers  as  being  in  excess  of  10^,  so  for  conical  bodies  the  transition 
region  will  begin  in  the  wake  and  laminar  flow  may  be  assumed.  At  lower 
altitudes,  say  around  100  kft,  the  transition  region  from  laminar  to  turbulent 
flow  would  begin  around  one  foot  and,  depending  on  the  type  of  vehicle,  it  would 
be  necessary  to  consider  turbulent  flow  over  the  body  at  these  altitudes.  This 
potential  development  has  been  kept  in  mind,  but  for  the  present,  only  laminar 
flow  is  considered. 

The  radiation  source  term  in  the  energy  equation  (3)  has  to  be  obtained  from  a 
radiation  transport  model.  Since  the  primary  focus  of  the  Phase  1  effort  has 
been  the  chemically  reacting  phenomena,  the  radiation  term  has  not  been  included 
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yet.  Under  the  proposed  Phase  II  effort,  radiation  models  would  be  reviewed  and 
an  appropriate  one  would  be  chosen  for  incorporation  into  the  code. 

In  summary,  the  system  of  equations  used  in  the  present  effort  consists  of 
the  partial  differential  equations  (1)  -  (3)  and  the  auxiliary  relations  (4)  - 
(11).  Given  the  mole  fractions  of  each  of  the  species  of  the  mixture  (or, 
equivalently,  given  the  mass  fractions),  these  equations  can  be  solved  for  the 
mixture  properties  p,  U,  h,  p,  and  T. 

The  numerical  method  used  to  solve  the  system  of  equations  is  based  on  the 
finite  difference  approach  discussed  by  Briley  and  McDonald  (Refs.  16,  17),  which 
uses  a  consistently  split,  linearized  block- implicit  (LBI)  ADI  procedure.  The  LBI 
scheme,  being  an  implicit  method,  does  not  suffer  from  the  stability  restrictions 
which  relate  the  temporal  step  to  the  spatial  mesh  size.  This  is  an  important 
advantage  in  view  of  the  existence  of  high  characteristic  velocities,  and  the  need 
to  use  locally  refined  meshes  for  accurate  prediction  of  the  flow  field  near  solid 
surfaces,  and  in  the  regions  of  sharp  temperature  or  species  concentration 
gradients.  Further,  the  solution  procedure  treats  the  nonlinearities 
noniterativelv  by  Taylor  series  linearization  to  the  requisite  order  in  time,  and 
then  splitting  the  matrix  into  a  series  of  easily  solved  block-banded  subsystems. 
The  solution  procedure  is,  thus,  computationally  efficient.  Further  details  on 
the  LBI  algorithm  are  available  in  Refs.  16  and  17. 

During  the  Phase  I  effort,  the  existing  scheme  for  the  solution  of  equations 
(1)  -  (3)  for  a  single  species  has  been  augmented  by  the  relations  (4)  -  (11)  to 
allow  the  solution  of  equations  (1)  -  (3)  for  a  gas  mixture,  given  the  mass 
fractions  (or  mole  fractions)  of  each  of  the  species  in  the  mixture. 


3.2.  Step  2  -  The  Species  Transport  Phase 

In  the  species  transport  phase,  each  species  is  represented  in  each  cell  by  a 
number  of  computational  particles.  A  computational  particle  can  be  thought  of  as 
a  collection  of  molecules  of  the  same  species  that  have  the  same  velocity.  Each 
computational  particle  is  then  moved  over  a  time  step  using  a  convection  velocity 
computed  in  the  continuum  phase,  and  an  appropriate  diffusion  velocity.  In  the 
DSMC  method,  this  diffusion  velocity  is  obtained  using  a  molecular  velocity 
distribution  function  (such  as  a  Maxwellian  distribution  function),  and  collisions 
between  the  "groups  of  molecules"  are  carried  out  at  the  end  of  each  time  step. 


In  the  present  approach,  a  macroscopic  diffusion  velocity  is  used,  rather  than  a 
molecular  one,  and  no  collisions  take  place  between  the  computational  particles. 
The  macroscopic  properties  that  result  from  the  collision  procedure  in  the  DSMC 
method,  such  as  the  viscosity  and  the  species  diffusivity  into  the  mixture,  are 
computed  explicitly,  and  the  motion  of  the  computational  particles  is  used  to 
simulate  the  solution  of  the  species  convection-diffusion  equations  by  means  of  a 
probabilistic  method.  The  advantage  of  the  present  approach  is  that  a  larger  time 
step  can  be  used  than  the  one  required  in  the  DSMC  method. 

The  two  aspects  of  the  species  transport  phase  that  need  to  be  discussed  are 
(a)  the  computation  of  the  macroscopic  diffusion  velocity  and  (b)  the  numerical 
scheme  for  moving  the  computational  particles . 


The  macroscopic  diffusion  velocity  of  a  particle  is  obtained  by  randomly 
picking  velocity  components  from  a  Gaussian  distribution  f(^)  with  zero  mean  and 
standard  deviation  (2DAt)^, 


f(a 


1  •-^^/(4DAt) 

2  (TtDAt)  ^  ® 


(12) 


where  At  is  the  time  step  over  which  the  particle  has  to  be  moved,  and  where  D  is 
the  diffusivity  of  the  species  under  consideration  into  the  mixture.  This  result 
follows  from  the  requirement  that  the  motion  of  a  large  number  of  particles  leads 
to  a  particle  distribution  (i.e.  species  concentration  distribution)  that “is 
consistent  with  the  solution  of  the  diffusion  equation  for  the  species  under 
consideration  (see,  for  example.  Refs.  18  and  19).  The  diffusivity  of  the 
species  i  into  the  mixture  can  be  determined  from  the  relation  (Ref.  13) 


z 

j=l 


(13) 


where  is  the  binary  diffusion  coefficient  of  species  i  into  species  j  ,  which 

can  be  expressed  as 
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(14) 


^ij  = 


r  1  11 

0.001858  t3/2 

^  J 

poij  0, 


Here  D^j  is  in  cm^/s,  p  in  atm,  T  in  K,  and  ^^iji  the  effective  collision  diameter 
given  by  =  h^o  +  Oj),  where  is  the  molecular  collision  diameter  of 

species  i,  in  A.  The  diffusion  collision  integral  fljj  can  be  approximated  by  (Ref. 
14). 


where 


T 


1 -0.145 


+ 


0.5 


'1-2.0 


(15) 


(T  being  the  molecular  temperature  diameter  of  species  i) .  The  values 
^  i 

of  a  and  T^  used  in  the  computation  of  the  dif fusivities  via  Eqs.  (14)  -  (16)  are 
the  same  as  those  used  in  the  computation  of  the  viscosities  via  Eqs.  (8)  -  (9). 
Once  the  proper  diffusivity  D  has  been  computed  using  Eqs.  (13)  -  (16),  diffusion 
velocity  components  can  be  picked  from  the  distribution  (12),  and  the  particle 
velocity  is  obtained  by  adding  the  diffusion  velocity  to  the  convection  velocity 
as  computed  in  the  continuum  phase . 


The  numerical  scheme  used  to  move  the  computational  particles  is  based  on  the 
particle  transport  model  incorporated  into  the  CELMINT  (Combined  Eulerian 
l^grangian  Multidimensional  Implicit  Navier-Stokes  Time -Dependent)  code  developed 
for  the  computation  of  flows  that  contain  solid  particles  or  liquid  droplets  (Ref. 
20) .  The  key  feature  of  this  particle  transport  model  is  that  it  integrates  the 
Lagrangian  equations  of  motion  for  a  particle  in  computational  space  rather  than 
physical  space.  This  eliminates  the  need  to  search  for  the  mesh  cell  in  which  a 
particle  resides.  The  present  application  of  this  model  differs  from  the  previous 
applications  of  the  SRA  code  (Refs.  21,  22)  in  that  there  is  no  force  acting  on 
the  computational  particle.  The  equation  of  motion  is,  in  this  case. 
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where  x  is  the  position  vector  of  the  particle,  and  V  is  the  particle  velocity. 
In  computational  space,  this  equation  becomes 


=  y  ^  V 

dt  axj  3 


(i  =  1,2,3) 


(18) 


where  are  the  computational  coordinates,  are  the  physical  coordinates  (the 
components  of  x) ,  and  are  the  physical  velocity  components  (the  components  of 
V)  of  the  particle. 

Integration  of  Eq.  (18)  allows  tracking  of  the  particle  motion  in  the 
computational  space.  As  mentioned  above,  this  eliminates  the  need  to  search  for 
the  mesh  cell  in  which  a  particle  resides.  In  addition,  this  simplifies  the 
search  for  boundaries,  needed  to  determine  whether  a  particle  has  reached  a  solid 
boundary  and  should  be  reflected,  or  a  particle  has  left  the  computational  domain 
altogether  and  should  be  omitted.  As  a  consequence,  the  method  can  be  applied 
easily  to  particle  motion  in  rather  complex  geometries  (Ref.  22) . 


In  summary,  the  species  transport  phase  consists  of  the  following  steps: 

(i)  Given  the  mass  fractions  (or  mole  fractions)  of  each  of  the  species  in  a 
mesh  cell,  the  number  of  computational  particles  of  each  species  in  each 
cell  is  determined. 

(ii)  The  computational  particles  are  distributed  randomly  across  the  mesh  cell. 

(iii)  Each  particle  is  moved  over  a  time  step  using  a  mixture  convection  velocity 
computed  in  the  continu\im  phase,  and  a  diffusion  velocity  randomly  picked 
from  a  distribution  based  on  the  species  dispersion  into  the  mixture  (Eqs. 
12-16)  . 

At  the  end  of  the  species  transport  phase,  a  new  distribution  of  computational 

particles  is  obtained,  to  be  used  in  the  reaction  phase,  described  below. 


3.3  Step  3  -  The  Chemical  Reaction  Phase 

In  the  reaction  phase,  a  representative  set  of  reactions  that  takes  place  in 
the  time  interval  At  is  computed  among  computational  particles  by  means  of 
techniques  similar  to  those  used  in  DSMC  methods.  Since  the  choice  of  molecular 
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collision/reaction  models  is  an  important  aspect  of  the  DSMC  methods,  a  brief 
history  is  presented  first. 

Early  applications  of  the  direct  simulation  methods  employed  the  well 
established  models  of  classical  kinetic  theory  (Refs.  12,  23).  The  hard  sphere 
model  is  the  easiest  to  apply,  and  both  it  and  the  Maxwell  model  were  widely  used 
for  early  studies  in  which  the  primary  objective  was  to  make  comparisons  with 
other  numerical  or  analytical  methods  employing  the  same  model.  These  models  are 
associated  with  unrealistic  viscosity- temperature  power  laws  and  became  inadequate 
when  serious  comparisons  were  to  be  made  with  experiment.  The  structure  of  strong 
shock  waves  was  the  subject  of  the  first  serious  comparisons  of  DSMC  results  and 
experiment.  The  inverse  power  law  of  repulsion  model  was  first  implemented  (Ref. 
24)  and  some  calculations  were  made  (Ref.  25)  with  the  Exp'^  and  Lennard- Jones 
models  that  include  the  long  range  attractive  force.  The  simpler  model  was  found 
to  be  generally  adequate  for  engineering  studies.  In  all  the  successful 
correlations  between  experiment  and  simulation,  the  molecular  impact  parameters 
were  chosen  to  match  the  coefficient  of  viscosity  of  the  real  gas. 

Experience  was  gradually  accumulated  on  the  effects  of  the  molecular  model  on 
representative  flows.  It  was  found  that  the  observed  effects  could  be  fully 
correlated  with  the  variation,  due  to  their  relative  velocity  in  the  collision,  of 
the  differential  cross-section  of  the  molecules.  The  form  of  the  deflection  angle 
scattering  law  of  the  molecules  appeared  to  be  comparatively  unimportant.  These 
considerations  led  to  the  introduction  (Ref.  26)  of  the  variable  hard  sphere,  or 
VHS ,  molecular  model .  This  is  essentially  a  hard  sphere  with  a  diameter  that 
varies  as  some  inverse  power  of  the  relative  velocity  in  the  collision.  The  VHS 
model  is  the  simplest  one  that  is  capable  of  modeling  the  viscosity  coefficient, 
including  its  temperature  dependence,  of  real  molecules.  It  has  been  found  (Refs. 
27,  18)  that,  as  long  as  this  is  done,  the  results  are  independent  of  the 
molecular  model. 

Nonequilibriura  chemical  reactions  are  readily  handled  by  procedures  that  are 
essentially  extensions  of  the  elementary  collision  theory  of  chemical  physics.  In 
this,  the  binary  reaction  rate  is  obtained  as  the  product  of  the  collision  rate 
for  collisions  with  energy  in  excess  of  the  activation  energy  and  the  probability 
of  reaction  or  steric  factor.  Since  the  formulation  is  in  terms  of  kinetic 
theory,  the  model  can  be  incorporated  directly  into  the  DSMC  method. 

The  chemical  data  for  gas  phase  reactions  is  almost  invariably  quoted  in  terms 
of  the  reaction  rate  constant.  A  form  of  the  collision  theory  that  is  consistent 
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with  the  VHS  model  can  be  used  (Ref.  7)  to  convert  the  temperature  dependent  rate 
constants  of  continuum  theory  into  collision  energy  dependent  steric  factors  of 
reaction  cross-sections.  Alternatively,  the  reactive  cross-sections  may  be 
available  directly  from  the  thermophysical  database  (Refs.  1,  5,  6,  25).  While 
this  is  more  in  keeping  with  the  simulation  approach,  it  leads  to  increased 
difficulties  in  ensuring  that  the  data  for  the  forward  and  reverse  reactions  is 
consistent  with  chemical  equilibrium  in  situations  where  the  reactions  should 
proceed  to  equilibrium.  Extensive  calculations  have  been  made  (Ref.  4)  for  the 
atmospheric  reentry  of  the  Space  Shuttle  and  direct  comparisons  were  made  with 
flight  data  and  with  overlapping  continurun  calculations.  Recombination  reactions 
are  not  important  at  very  low  densities  where  the  probability  of  three -body 
collisions  is  extremely  small.  They  become  important  in  the  continuum  overlap 
region  and  these  are  dealt  with  by  assigning  a  lifetime  to  each  binary  collision 
and  regarding  the  ternary  collision  as  a  further  binary  collision  between  the  pair 
of  molecules  in  the  original  collision  and  a  third  molecule.  These  procedures  can 
be  incorporated  into  the  DSMC  method  in  a  form  (Ref.  29)  that  guarantees  the 
correct  equilibrium  state. 

The  present  approach  follows  the  procedure  described  by  Bird  (Ref.  7),  which 
uses  the  VHS  molecular  model  and  reaction  cross-sections  that  are  based  on  the 
(temperature  dependent)  rate  constants  of  continuum  theory.  The  numerical  scheme 
consists  of  the  following  steps  (for  each  mesh  cell): 

(i)  Pick  a  reaction  among  the  set  of  all  binary  reactions  under  consideration, 
say  A  +  B  -»  C  +  D. 

(ii)  Pick  one  particle  of  each  of  the  (two)  species  A  and  B.  Note  that  the 
statistical  weight  of  these  particles  (i.e.  the  number  of  physical 
molecules  that  they  represent)  must  be  the  same.  In  the  DSMC  method,  the 
particles  chosen  have  velocities  associated  with  them.  In  the  present 
method,  the  choice  of  a  particle  includes  the  choice  of  velocity 
components  from  the  appropriate  (for  example,  Maxwellian)  molecular 
velocity  distribution. 

(iii)  Determine  the  reaction  cross-section. 
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(iv)  Based  on  the  reaction  probability,  which  relates  directly  to  the  reaction 
cross-section,  determine  whether  or  not  a  reaction  takes  place. 

(v)  If  a  reaction  takes  place,  create  particles  of  the  species  C  and  D,  and 
omit  the  particles  of  species  A  and  B.  In  other  words,  update  the  species 
concentrations.  In  addition,  update  a  reaction  time  counter  by  an  amount 
inversely  proportional  to  the  reaction  probability. 

(vi)  Omit  the  reaction  from  the  set  of  reactions  under  consideration  if  its 
time  counter  exceeds  the  physical  time. 

(vii)  Repeat  the  procedure  until  the  set  of  reactions  is  empty,  i.e.  until  all 
reaction  time  counters  have  exceeded  the  physical  time. 

The  use  of  reaction  time  counters  ensures  that,  over  a  large  number  of 
reactions,  the  procedure  leads  to  a  reaction  rate  that  is  in  agreement  with  the 
reaction  cross-section  (cf.  Ref.  7).  This  leads  to  an  interesting  observation: 
if  the  reaction  cross-sections  are  determined  such  that  the  resulting  reaction 
rates  correspond  to  those  of  continuum  theory,  then  the  steps  (ii)  -  (v)  in  the 
above  procedure  could  be  replaced  by 

(v)  Let  a  reaction  take  place.  Update  the  species  concentrations;  update  the 
reaction  time  counter  by  an  amount  inversely  proportional  to  the  reaction 
rate . 


In  that  case,  the  procedure  can  be  considered  as  a  stochastic  method  for 
solving  a  system  of  ordinary  differential  equations  of  the  form 


dt 


B  C,D 


(19) 


where  n^  denotes  the  number  density  of  species  A,  and  k^g(T)  and  kQf)(T)  are  the 
(continuum)  rate  coefficients.  The  first  sum  on  the  right-hand  side  of  Eq.  (19) 
is  over  all  reaction  of  the  form  A  +  B  -*  . . . ,  while  the  second  sum  is  over  all 
reactions  of  the  form  C  +  D  -*  A  +  ... 


-16- 


The  advantage  of  using  the  complete  procedure  (i)  -  (vii)  rather  than  the 
simplified  one  (with  (ii)  -  (v)  replaced  by  (v)*)  lies  in  the  fact  that 
alternative  formulations  for  the  reaction  cross-sections,  based  on  molecular 
models  rather  than  continuum  theory,  can  be  implemented  easily. 


3.4.  Coupling  of  the  Three  Phases 

Finally,  it  is  necessary  to  make  some  remarks  about  the  coupling  of  the  three 
phases  in  the  hybrid  scheme  developed  under  the  present  effort.  The  continuum 
phase  provides  the  species  transport  phase  with  a  convection  velocity  and  with 
mixture  properties  that  are  used  to  compute  the  diffusion  velocity.  The  species 
mass  fractions  (or  mole  fractions),  which  are  not  changed  by  the  continuum  phase, 
and  the  mixture  density  provide  the  species  number  densities,  used  to  determine 
the  number  of  computational  particles  of  each  species  and  their  statistical 
weights . 

The  species  transport  phase  moves  the  computational  particles ,  thereby 
updating  the  species  number  densities  in  each  mesh  cell.  The  chemical  reaction 
phase  updates  the  number  densities  once  more  when  reactions  take  place.  This 
phase  also  uses  some  continuum  properties  (such  as  temperature)  in  the  computation 
of  reaction  cross-sections.  At  the  end  of  the  reaction  phase,  new  species  mass 
fractions  (or  mole  fractions)  can  be  computed  in  each  mesh  cell.  As  such,  the 
combined  transport/reaction  phase  serves  to  update  only  the  species  mass 
fractions . 

A  word  of  caution  is  necessary  here.  The  species  mass  fractions  that  are 
computed  at  the  end  of  the  combined  transport/reaction  phase  make  sense  in  the 
continuum  phase  only  if  a  very  large  number  of  computational  particles  is  used, 
or,  equivalently,  if  an  ensemble  average  is  taken  over  a  large  number  of  separate 
calculations.  Since  for  steady-state  calculations  a  time-average  is  equivalent  to 
an  ensemble  average,  it  is  advantageous  to  reduce  the  number  of  particles  used  in 
each  time  step,  and  to  time -average  the  mass  fractions  computed  by  the 
transport/reaction  phase  before  they  are  used  in  the  continuum  phase.  Due  to  the 
coupling  of  the  continuum  phase  and  the  transport/reaction  phase  this  cannot  be 
done  efficiently  directly,  but  as  can  be  shown,  the  required  averaging  can  be 
accomplished  via  underrelaxation  of  the  changes  in  the  species  mass  fractions 
before  they  are  used  in  the  continuum  phase. 
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4.  RESULTS  AND  CONCLUSIONS 


The  hybrid  approach  described  above  has  been  developed  under  the  Phase  I 
program.  In  order  to  demonstrate  the  feasibility  of  this  approach,  a  calculation 
has  been  performed  for  the  axisymmetric  hypersonic  laminar  flow  of  dissociating 
air  over  a  10  -  degree  cone.  Although  this  represents  a  relatively  simple  test 
case  chosen  to  clarify  the  evaluation  of  the  procedure,  it  does  contain  all  the 
basic  elements  of  hypersonic  vehicle  environment  simulation  insofar  as  they  have 
been  included  in  the  present  model. 

The  relevant  flow  conditions  for  the  test  case  are: 

Altitude  45.72  km  (150  kft) 

Velocity  3.27  km/s  (10.73  kft/s) 

Mach  Number  10 

Ambient  Pressure  136.7  N/m^ 

Ambient  Temperature  266.2  K 

Ambient  Mole  Fractions: 

O2  0.24 

N2  0.76 

The  wall  was  assumed  adiabatic  and  non-catalytic .  The  chemical  reactions 
considered  in  the  present  study  consist  of  19  reactions  for  5  species  (N2,  O2 ,  NO, 
N  and  0),  and  the  data  are  listed  in  Table  1,  which  was  taken  from  Ref.  30. 
Three-body  atom  recombination  reactions  are  not  significant  at  high  altitudes,  and 
have  not  been  included. 

Figure  1  shows  the  computational  domain  and  the  numerical  grid  used  in  the 
calculations.  The  downstream  boundary  of  the  computational  domain  (boundary  DE  in 
Fig.  1)  is  about  3  inches  from  the  vertex  of  the  cone  (point  C) ,  to  ensure  laminar 
flow  over  the  section  CD  of  the  cone  under  consideration.  The  upstream  boundary 
AB  and  the  "far -field"  boundary  AE  were  chosen  such  that  they  were  both  inflow 
boundaries.  The  bow  shock  wave  intersects  the  downstream  boundary.  The  grid  used 
contains  100  x  150  mesh  points,  concentrated  near  the  cone  surface  and  the  cone 
vertex,  to  allow  proper  resolution  of  the  boundary  layer. 

Results  of  the  computation  are  shown  in  Figures  2-8.  Figure  2  shows  contours 
of  the  density,  the  pressure,  and  the  temperature.  In  this  figure,  the  shock  wave 
and  the  bou  dary  layer  can  be  distinguished  easily.  The  location  and  strength  of 
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the  bow  shock  are  in  good  agreement  with  the  inviscid  theory  of  conical  flow,  even 
though  no  particular  attention  was  paid  to  resolving  the  bow  shock.  As  can  be 
seen  in  Figure  2b,  the  pressure  is  approximately  constant  across  the  boundary 
layer.  In  Figures  3-5,  the  temperature  calculated  for  the  reacting  flow  is 
compared  with  the  temperature  that  can  be  obtained  for  non- reacting  flow,  without 
species  transport,  i.e.  for  the  flow  with  "frozen"  species  concentrations.  Figure 
3  shows  that  the  temperature  on  the  cone  surface  is  lower  for  the  reacting  flow 
than  for  the  non- reacting  flow.  This  effect  is  to  be  expected,  because  the 
formation  of  0,  NO,  and  N  requires  energy.  The  high  temperature  calculated  near 
the  cone  vertex  indicates  that  the  flow  near  this  vertex  has  been  resolved 
adequately.  In  Figure  4,  the  temperature  profile  is  shown  at  an  axial  distance  of 
2.5  inches  downstream  of  the  cone  vertex.  The  thickness  of  the  shock  indicates  a 
lack  of  grid  resolution  in  the  shock  wave.  Again,  the  temperature  of  the  reacting 
flow  is  seen  to  be  lower  than  the  temperature  of  the  non-reacting  flow.  Figure  5a 
is  an  enlargement  of  Figure  4  in  the  boundary  layer  region,  while  Figures  5b-c 
show  the  temperature  profiles  at  axial  distances  of  1.0  and  0.1  inch  downstream  of 
the  cone  vertex. 

Figure  6  shows  contours  of  the  mass  fractions  of  0,  NO,  and  N  in  the  region 
from  2.0  to  2.5  inches  downstream  of  the  cone  vertex.  The  mass  fractions  can  be 
seen  to  increase  in  the  axial  direction  because  0,  NO  and  N  are  not  only  produced 
by  the  chemical  reactions,  but  also  convected  downstream  by  the  mean  flow.  The 
increase  of  the  mass  fractions  along  the  cone  surface  is  also  clearly  visible  in 
Figure  7.  Finally,  Figure  8  shows  profiles  of  the  mass  fractions  at  several  axial 
locations . 

The  results  of  the  calculation  discussed  above  demonstrate  that  the  hybrid 
scheme  is  capable  of  computing  chemically  reacting  flows  and  has  the  potential 
for  efficiently  computing  high  Mach  number  flows  around  vehicles  in  the  upper 
layers  of  the  atmosphere.  Further  development  of  this  scheme  is  the  focus  of  the 
proposed  Phase  II  effort. 
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Table  1.  List  of  Chemical  Reactions. 

Five -species  Model.® 

Rate  Coefficient,  in  m^/molecule— s 
aT^  exp(-c/T) 

Number  Reaction  a  be 


1 

^2 

+ 

N 

20 

+ 

N 

5.993x10"12 

0 

• 

1 

59370 

2 

^2 

+ 

NO 

- 

20 

+ 

NO 

5.993x10“12 

-1.0 

59370 

3 

02 

+ 

N2 

- 

20 

+ 

N2 

1.198xl0“ll 

-1.0 

59370 

4 

02 

+ 

O2 

- 

20 

+ 

O2 

5.393x10“^! 

-1.0 

59370 

5 

O2 

+ 

0 

20 

+ 

0 

1.498x10“^° 

-1.0 

59370 

6 

N2 

+ 

0 

2N 

+ 

0 

3.187x10“13 

-0.5 

113000 

7 

N2 

+ 

02 

2N 

+ 

02 

3.187x10“13 

-0.5 

113000 

8 

N2 

+ 

NO 

- 

2N 

+ 

NO 

3.187x10“13 

-0.5 

113000 

9 

N2 

-t- 

N2 

- 

2N 

+ 

N2 

7.968x10“^3 

-0.5 

113000 

10 

N2 

+ 

N 

2N 

+ 

N 

6.900x10"® 

-1.5 

113000 

11 

NO 

+ 

N2 

N 

+ 

0 

-F 

N2 

6.590xl0"10 

-1.5 

75550 

12 

NO 

+ 

O2 

-♦ 

N 

+ 

0 

+ 

O2 

6.590xl0"10 

-1.5 

75550 

13 

NO 

+ 

NO 

N 

+ 

0 

+ 

NO 

1.318x10"® 

-1.5 

75550 

14 

NO 

+ 

0 

N 

+ 

0 

+ 

0 

1.318x10"® 

-1.5 

75550 

15 

NO 

+ 

N 

N 

+ 

0 

+ 

N 

1. 318x10"® 

-1.5 

75550 

16 

NO 

+ 

0 

- 

O2 

+ 

N 

5.279x10"21 

1.0 

19700 

17 

N2 

+ 

0 

NO 

+ 

N 

1. 120x10"!® 

0.0 

37500 

18 

O2 

+ 

N 

- 

NO 

+ 

0 

1.598x10"!® 

0.5 

3600 

19 

NO 

+ 

N 

- 

N2 

+ 

0 

2.490x10"!'^ 

0.0 

0 

®  Without  three-body  recombination  reactions  which  do  not  participate 
at  low  densities. 
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;ure  1.  Computational  Domain  an 


(a)  Density 


(b)  Pressure 


(c)  Temperature 


Figure  2.  Contours  of  Density,  Pressure,  and  Temperature 


Between  x 


0  a  n  d  X  =  1 
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Figure  4,  Radial  Variation  of  the' Temperature  at 


-2  5- 


5a,  Radial  Variation  of  the  Temperature  at 


Radial  Variation  of  the  Temperature 


Radial 


Mass  Fraction 
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Figure  7.  Axial  Variation  of  the  Mass  Fractions  Along  the  Cone  Surface 
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Figure  8a.  Radial  Variation  of  the  Mass  Fractio 


Radial  Variation  of  the  Mass  Frac 


5 .  FUTURE  DEVELOPMENT 


As  described  above,  under  the  Phase  I  effort  a  hybrid  Navier-Stokes/Monte 
Carlo  method  has  been  developed  for  the  computation  of  hypersonic,  viscous, 
chemically  reacting  flows.  This  method  consists  of  three  major  phases.  The  first 
phase  is  the  continuum  phase  where  the  mean  flow  field  calculations  are  performed 
via  numerical  solution  of  the  Navier-Stokes  equations  and  the  global  continuity 
and  energy  equations  for  the  gas  mixture.  The  second  phase  is  the  species 
transport  phase,  where  the  "computational  particles"  representing  the  species 
molecules  are  transported  simulating  the  convection  and  diffusion  of  the  species. 
Finally,  the  third  phase  is  the  reaction  phase,  where  appropriate  collision  models 
are  used  in  a  Monte  Carlo-like  way  to  simulate  the  chemical  kinetics.  The  first 
phase  uses  an  Eulerian  description,  while  the  second  and  third  phases  utilize  a 
statistical  Lagrangian  description.  The  method  offers  the  flexibility  of  the 
DSMC-method  with  respect  to  the  chemical  reactions,  but  allows  the  use  of  larger 
time  steps,  because  molecular  collision  properties  such  as  viscosity  and 
diffusivity  are  obtained  from  semi-empirical  relations. 

A  future  effort  would  concentrate  on  the  further  development  of  both  the 
numerical  and  the  physical  aspects  of  the  hybrid  scheme,  with  the  objective  of 
developing  and  demonstrating  a  computer  code  that  can  be  used  as  a  tool  for 
investigating  vehicle  performance  and  sensor  visibility  in  hypersonic  flight 
within  the  atmosphere. 

Specific  tasks  that  would  comprise  this  effort  are  the  following: 

(i)  Further  development  of  the  Lagrangian  statistical  species  transport 
algorithm  via  the  inclusion  of  sophisticated  statistical  models. 

(ii)  Added  capabilities  in  terms  of  boundary  conditions  (catalytic  walls), 
recombination  reactions,  and  ionization. 

(iii)  Addition  of  a  radiation  transport  model. 

(iv)  Demonstration  of  the  fully  developed  procedure  via  one  or  more  relevant 
flow  field  calculations. 
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(v)  Development  of  suitable  input  and  output  packages  for  use  of  the  code  in 

a  workstation  environment. 

These  tasks  have  been  detailed  in  a  Phase  II  proposal. 
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